TODO
options(timeout = 10000)
dir.create("data/steinbock/raw", recursive = TRUE)
download.file("https://zenodo.org/record/6642699/files/panel.csv",
"data/steinbock/panel.csv")
download.file("https://zenodo.org/record/5949116/files/Patient1.zip",
"data/steinbock/raw/Patient1.zip")
download.file("https://zenodo.org/record/5949116/files/Patient2.zip",
"data/steinbock/raw/Patient2.zip")
download.file("https://zenodo.org/record/5949116/files/Patient3.zip",
"data/steinbock/raw/Patient3.zip")
download.file("https://zenodo.org/record/5949116/files/Patient4.zip",
"data/steinbock/raw/Patient4.zip")
download.file("https://zenodo.org/record/5949116/files/compensation.zip",
"data/compensation.zip")
unzip("data/compensation.zip", exdir="data", overwrite=TRUE)
unlink("data/compensation.zip")
download.file("https://zenodo.org/record/5949116/files/sample_metadata.xlsx",
destfile = "data/sample_metadata.xlsx")
download.file("https://zenodo.org/record/7079294/files/gated_cells.zip",
"data/gated_cells.zip")
unzip("data/gated_cells.zip", exdir="data", overwrite=TRUE)
unlink("data/gated_cells.zip")
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
devtools::install_github(c("BodenmillerGroup/imcRtools",
"BodenmillerGroup/cytomapper"))
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("pheatmap", "viridis",
"tiff", "distill", "openxlsx", "ggrepel", "patchwork",
"mclust", "RColorBrewer", "uwot", "Rtsne", "caret",
"randomForest", "ggridges", "gridGraphics", "scales",
"CATALYST", "scuttle", "scater", "dittoSeq",
"tidyverse", "batchelor", "bluster","scran"))
steinbock preprocess imc images --hpf 50
## bash: no job control in this shell
## 2022-09-25 17:03:57,429 INFO steinbock - img/Patient4_005.tiff
## 2022-09-25 17:03:58,802 INFO steinbock - img/Patient4_006.tiff
## 2022-09-25 17:04:00,339 INFO steinbock - img/Patient4_007.tiff
## 2022-09-25 17:04:02,644 INFO steinbock - img/Patient4_008.tiff
## 2022-09-25 17:04:03,881 INFO steinbock - img/Patient3_001.tiff
## 2022-09-25 17:04:05,098 INFO steinbock - img/Patient3_002.tiff
## 2022-09-25 17:04:06,099 INFO steinbock - img/Patient3_003.tiff
## 2022-09-25 17:04:07,314 INFO steinbock - img/Patient2_001.tiff
## 2022-09-25 17:04:08,244 INFO steinbock - img/Patient2_002.tiff
## 2022-09-25 17:04:09,352 INFO steinbock - img/Patient2_003.tiff
## 2022-09-25 17:04:10,534 INFO steinbock - img/Patient2_004.tiff
## 2022-09-25 17:04:12,851 INFO steinbock - img/Patient1_001.tiff
## 2022-09-25 17:04:14,181 INFO steinbock - img/Patient1_002.tiff
## 2022-09-25 17:04:15,607 INFO steinbock - img/Patient1_003.tiff
## 2022-09-25 17:04:32,808 INFO steinbock - images.csv
The step took 1.9 minutes.
steinbock segment deepcell --minmax
## bash: no job control in this shell
## /opt/venv/lib/python3.8/site-packages/deepcell_toolbox/deep_watershed.py:179: FutureWarning: `selem` is a deprecated argument name for `h_maxima`. It will be removed in version 1.0. Please use `footprint` instead.
## markers = h_maxima(image=maxima,
## 2022-09-25 17:05:28,489 INFO steinbock - masks/Patient1_001.tiff
## 2022-09-25 17:05:45,261 INFO steinbock - masks/Patient1_002.tiff
## 2022-09-25 17:06:02,687 INFO steinbock - masks/Patient1_003.tiff
## 2022-09-25 17:06:19,024 INFO steinbock - masks/Patient2_001.tiff
## 2022-09-25 17:06:37,602 INFO steinbock - masks/Patient2_002.tiff
## 2022-09-25 17:07:07,558 INFO steinbock - masks/Patient2_003.tiff
## 2022-09-25 17:07:26,997 INFO steinbock - masks/Patient2_004.tiff
## 2022-09-25 17:07:44,841 INFO steinbock - masks/Patient3_001.tiff
## 2022-09-25 17:08:07,048 INFO steinbock - masks/Patient3_002.tiff
## 2022-09-25 17:08:23,965 INFO steinbock - masks/Patient3_003.tiff
## 2022-09-25 17:08:38,329 INFO steinbock - masks/Patient4_005.tiff
## 2022-09-25 17:08:55,132 INFO steinbock - masks/Patient4_006.tiff
## 2022-09-25 17:09:12,511 INFO steinbock - masks/Patient4_007.tiff
## 2022-09-25 17:09:27,149 INFO steinbock - masks/Patient4_008.tiff
The step took 4.91 minutes.
steinbock measure intensities
## bash: no job control in this shell
## 2022-09-25 17:09:35,277 INFO steinbock - intensities/Patient1_001.csv
## 2022-09-25 17:09:36,793 INFO steinbock - intensities/Patient1_002.csv
## 2022-09-25 17:09:38,147 INFO steinbock - intensities/Patient1_003.csv
## 2022-09-25 17:09:39,407 INFO steinbock - intensities/Patient2_001.csv
## 2022-09-25 17:09:41,371 INFO steinbock - intensities/Patient2_002.csv
## 2022-09-25 17:09:42,629 INFO steinbock - intensities/Patient2_003.csv
## 2022-09-25 17:09:44,325 INFO steinbock - intensities/Patient2_004.csv
## 2022-09-25 17:09:46,051 INFO steinbock - intensities/Patient3_001.csv
## 2022-09-25 17:09:47,694 INFO steinbock - intensities/Patient3_002.csv
## 2022-09-25 17:09:49,691 INFO steinbock - intensities/Patient3_003.csv
## 2022-09-25 17:09:52,039 INFO steinbock - intensities/Patient4_005.csv
## 2022-09-25 17:09:55,301 INFO steinbock - intensities/Patient4_006.csv
## 2022-09-25 17:09:58,019 INFO steinbock - intensities/Patient4_007.csv
## 2022-09-25 17:10:00,812 INFO steinbock - intensities/Patient4_008.csv
The step took 0.49 minutes.
steinbock measure regionprops
## bash: no job control in this shell
## 2022-09-25 17:10:04,682 INFO steinbock - regionprops/Patient1_001.csv
## 2022-09-25 17:10:06,010 INFO steinbock - regionprops/Patient1_002.csv
## 2022-09-25 17:10:07,476 INFO steinbock - regionprops/Patient1_003.csv
## 2022-09-25 17:10:10,627 INFO steinbock - regionprops/Patient2_001.csv
## 2022-09-25 17:10:13,100 INFO steinbock - regionprops/Patient2_002.csv
## 2022-09-25 17:10:14,512 INFO steinbock - regionprops/Patient2_003.csv
## 2022-09-25 17:10:16,217 INFO steinbock - regionprops/Patient2_004.csv
## 2022-09-25 17:10:17,885 INFO steinbock - regionprops/Patient3_001.csv
## 2022-09-25 17:10:19,749 INFO steinbock - regionprops/Patient3_002.csv
## 2022-09-25 17:10:21,819 INFO steinbock - regionprops/Patient3_003.csv
## 2022-09-25 17:10:23,363 INFO steinbock - regionprops/Patient4_005.csv
## 2022-09-25 17:10:25,371 INFO steinbock - regionprops/Patient4_006.csv
## 2022-09-25 17:10:26,722 INFO steinbock - regionprops/Patient4_007.csv
## 2022-09-25 17:10:27,960 INFO steinbock - regionprops/Patient4_008.csv
The step took 0.45 minutes.
steinbock measure neighbors --type expansion --dmax 4
## bash: no job control in this shell
## 2022-09-25 17:10:33,154 INFO steinbock - neighbors/Patient1_001.csv
## 2022-09-25 17:10:36,080 INFO steinbock - neighbors/Patient1_002.csv
## 2022-09-25 17:10:39,505 INFO steinbock - neighbors/Patient1_003.csv
## 2022-09-25 17:10:42,409 INFO steinbock - neighbors/Patient2_001.csv
## 2022-09-25 17:10:45,442 INFO steinbock - neighbors/Patient2_002.csv
## 2022-09-25 17:10:47,932 INFO steinbock - neighbors/Patient2_003.csv
## 2022-09-25 17:10:51,513 INFO steinbock - neighbors/Patient2_004.csv
## 2022-09-25 17:10:55,454 INFO steinbock - neighbors/Patient3_001.csv
## 2022-09-25 17:10:58,535 INFO steinbock - neighbors/Patient3_002.csv
## 2022-09-25 17:11:02,234 INFO steinbock - neighbors/Patient3_003.csv
## 2022-09-25 17:11:05,598 INFO steinbock - neighbors/Patient4_005.csv
## 2022-09-25 17:11:11,984 INFO steinbock - neighbors/Patient4_006.csv
## 2022-09-25 17:11:15,531 INFO steinbock - neighbors/Patient4_007.csv
## 2022-09-25 17:11:19,859 INFO steinbock - neighbors/Patient4_008.csv
The step took 0.87 minutes.
library(imcRtools)
spe <- read_steinbock("data/steinbock/")
The step took 0.56 minutes.
library(openxlsx)
library(tidyverse)
colnames(spe) <- paste0(spe$sample_id, "_", spe$ObjectNumber)
meta <- read.xlsx("data/sample_metadata.xlsx")
spe$patient_id <- as.vector(str_extract_all(spe$sample_id, "Patient[1-4]",
simplify = TRUE))
spe$ROI <- as.vector(str_extract_all(spe$sample_id, "00[1-8]",
simplify = TRUE))
spe$indication <- meta$Indication[match(spe$patient_id, meta$Sample.ID)]
rowData(spe)$use_channel <- !grepl("DNA|Histone", rownames(spe))
assay(spe, "exprs") <- asinh(counts(spe)/1)
The step took 0.02 minutes.
library(cytomapper)
images <- loadImages("data/steinbock/img/")
channelNames(images) <- rownames(spe)
The step took 0.24 minutes.
masks <- loadImages("data/steinbock/masks/", as.is = TRUE)
## All files in the provided location will be read in.
The step took 0 minutes.
patient_id <- str_extract_all(names(images), "Patient[1-4]", simplify = TRUE)
indication <- meta$Indication[match(patient_id, meta$Sample.ID)]
mcols(images) <- mcols(masks) <- DataFrame(sample_id = names(images),
patient_id = patient_id,
indication = indication)
The step took 0 minutes.
cytomapper_sce <- measureObjects(masks, image = images, img_id = "sample_id")
The step took 7.18 minutes.
sce <- readSCEfromTXT("data/compensation/")
## Spotted channels: Y89, In113, In115, Pr141, Nd142, Nd143, Nd144, Nd145, Nd146, Sm147, Nd148, Sm149, Nd150, Eu151, Sm152, Eu153, Sm154, Gd155, Gd156, Gd158, Tb159, Gd160, Dy161, Dy162, Dy163, Dy164, Ho165, Er166, Er167, Er168, Tm169, Er170, Yb171, Yb172, Yb173, Yb174, Lu175, Yb176
## Acquired channels: Ar80, Y89, In113, In115, Xe131, Xe134, Ba136, La138, Pr141, Nd142, Nd143, Nd144, Nd145, Nd146, Sm147, Nd148, Sm149, Nd150, Eu151, Sm152, Eu153, Sm154, Gd155, Gd156, Gd158, Tb159, Gd160, Dy161, Dy162, Dy163, Dy164, Ho165, Er166, Er167, Er168, Tm169, Er170, Yb171, Yb172, Yb173, Yb174, Lu175, Yb176, Ir191, Ir193, Pt196, Pb206
## Channels spotted but not acquired:
## Channels acquired but not spotted: Ar80, Xe131, Xe134, Ba136, La138, Ir191, Ir193, Pt196, Pb206
assay(sce, "exprs") <- asinh(counts(sce)/5)
The step took 0.06 minutes.
plotSpotHeatmap(sce)
sce2 <- binAcrossPixels(sce, bin_size = 10)
The step took 0.21 minutes.
library(CATALYST)
bc_key <- as.numeric(unique(sce$sample_mass))
bc_key <- bc_key[order(bc_key)]
sce <- assignPrelim(sce, bc_key = bc_key)
## Debarcoding data...
## o ordering
## o classifying events
## Normalizing...
## Computing deltas...
sce <- estCutoffs(sce)
sce <- applyCutoffs(sce)
sce <- filterPixels(sce, minevents = 40, correct_pixels = TRUE)
The step took 0.17 minutes.
sce <- computeSpillmat(sce)
sm <- metadata(sce)$spillover_matrix
The step took 0.03 minutes.
library(dittoSeq)
library(patchwork)
rowData(spe)$channel_name <- paste0(rowData(spe)$channel, "Di")
isotope_list <- CATALYST::isotope_list
isotope_list$Ar <- 80
spe <- compCytof(spe, sm,
transform = TRUE, cofactor = 1,
isotope_list = isotope_list,
overwrite = FALSE)
before <- dittoScatterPlot(spe, x.var = "Ecad", y.var = "CD303",
assay.x = "exprs", assay.y = "exprs") +
ggtitle("Before compensation")
after <- dittoScatterPlot(spe, x.var = "Ecad", y.var = "CD303",
assay.x = "compexprs", assay.y = "compexprs") +
ggtitle("After compensation")
before + after
assay(spe, "counts") <- assay(spe, "compcounts")
assay(spe, "exprs") <- assay(spe, "compexprs")
assay(spe, "compcounts") <- assay(spe, "compexprs") <- NULL
channelNames(images) <- rowData(spe)$channel_name
adapted_sm <- adaptSpillmat(sm, paste0(rowData(spe)$channel, "Di"),
isotope_list = isotope_list)
## Compensation is likely to be inaccurate.
## Spill values for the following interactions
## have not been estimated:
## Ir191Di -> Ir193Di
## Ir193Di -> Ir191Di
images_comp <- compImage(images, adapted_sm)
plotPixels(images[5], colour_by = "Yb173Di",
image_title = list(text = "Yb173 (Ecad) - before",
position = "topleft"),
legend = NULL, bcg = list(Yb173Di = c(0, 4, 1)))
plotPixels(images[5], colour_by = "Yb174Di",
image_title = list(text = "Yb174 (CD303) - before",
position = "topleft"),
legend = NULL, bcg = list(Yb174Di = c(0, 4, 1)))
plotPixels(images_comp[5], colour_by = "Yb173Di",
image_title = list(text = "Yb173 (Ecad) - after",
position = "topleft"),
legend = NULL, bcg = list(Yb173Di = c(0, 4, 1)))
plotPixels(images_comp[5], colour_by = "Yb174Di",
image_title = list(text = "Yb174 (CD303) - after",
position = "topleft"),
legend = NULL, bcg = list(Yb174Di = c(0, 4, 1)))
channelNames(images_comp) <- rownames(spe)
The step took 8.92 minutes.
set.seed(20220118)
img_ids <- sample(seq_len(length(images_comp)), 3)
cur_images <- images_comp[img_ids]
cur_images <- normalize(cur_images, separateImages = TRUE)
cur_images <- normalize(cur_images, inputRange = c(0, 0.2))
plotPixels(cur_images,
mask = masks[img_ids],
img_id = "sample_id",
missing_colour = "white",
colour_by = c("CD163", "CD20", "CD3", "Ecad", "DNA1"),
colour = list(CD163 = c("black", "yellow"),
CD20 = c("black", "red"),
CD3 = c("black", "green"),
Ecad = c("black", "cyan"),
DNA1 = c("black", "blue")),
image_title = NULL,
legend = list(colour_by.title.cex = 0.9,
colour_by.labels.cex = 0.9))